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ABSTRACT 

We present numerical models for supernova remnant evolution, using a new version of the hydrodynamical code SUPREMNA. We 
added cosmic ray diffusion equation to the code scheme, employing two-fluid approximation. We investigate the dynamics of the 
simulated supernova renmants with different values of cosmic ray acceleration efficiency and diffusion coefficient. We compare the 
numerical models with observational data of Tycho's and SN1006 supernova remnants. We find models which reproduce the observed 
locations of the blast wave, contact discontinuity, and reverse shock for the both remnants, thus allowing us to estimate the contribution 
of cosmic ray particles into total pressure and cosmic-ray energy losses in these supernova remnants. We derive that the energy losses 
due to cosmic rays escape in Tycho's supernova remnant are 10-20% of the kinetic energy flux and 20-50% in SN1006. 

Key words, acceleration of particles — diffusion — hydrodynamics — shock waves — Methods: numerical — ISM: supernova 
renmants 



1. Introduction 

Clear evidence for effective acceleration of cosmic-ray (CR) par- 
ticles in young supernova remnants (SNR) stems from TeV ob- 
servations of the Galactic sources by HESS (Aharonian et al. 
2005, 2006), CANGAROO-n (Katagiri et al. 2005), MAGIC 
(Albert et al. 2007), SHALON (Sinitsyna et al. 2009), VERITAS 
(Acciari et al. 2009, 2010). In addition, discoveries of non- 
thermal X-ray emission from SNRs (see a review by Reynolds 
2008) point to the presence of electrons accelerated to TeV en- 
ergies in supernova remnants (SNR). 

In the last decades a number of numerical methods which 
account for CR acceleration in simulations of supernova rem- 
nants were developed. An extensive review of some of these 
techniques can be found in Malkov & O'C Drury (2001) and 
Caprioli et al. (2010). Nearby Galactic SNRs provide an excel- 
lent opportunity to test these models and to study the efficiency 
of CR acceleration processes. For example, the proximity of the 
blast wave (BW) to contact discontinuity (CD) in Tycho SNR 
measured by Warren et al. (2005) is inconsistent with adiabatic 
hydrodynamic models of SNR evolution, and can be explained 
orily if cosmic-ray acceleration of the particles occurs at the for- 
ward shock. The similar evidence was presented for SN1006 
supernova remnant by Cassam-Chenai et al. (2008) and Miceli 
et al. (2009). 

Both these objects were already used in testing of the numer- 
ical and analytical models (e.g. Ellison 2001, Ellison et al. 2007, 
Volk et al. 2008) of CR acceleration in SNRs. 2D simulations of 
Tycho's SNR evolution and investigation of Rayleigh-Taylor in- 
stability development with the gas adiabatic index values down 
to y = 1.1 were performed by Wang (2010). 3D hydrodynamical 
modeling of Tycho was conducted by Ferrand et al. (2010). 



The effects of the shock modification by cosmic rays in 
Kepler SNR were studied by Decourchelle et al. (2000), who 
modeled the X-ray spectra using a non-linear non-equilibrium 
ionization method. A more detailed study of the acceleration ef- 
fects on the thermal emission from the shocked supernova ejecta 
was conducted by Patnaude et al. (2010). 

In this study we present hydrodynamical (HD) simulations 
of supernova remnant evolution with the account for diffusive 
cosmic-ray acceleration. We introduced a CR diffusion equation 
into the numerical code supremna, developed by Sorokina et al. 
(2004), Kosenko (2006). This code calculates the evolution of 
a supemova remnant assuming spherical symmetry and taking 
into account time-dependent ionization and thermal conduction. 
To include the effects of CR acceleration into the scheme, we ap- 
ply a two-fluid approximation, i.e. we introduce a CR diffusion 
equation into the system of hydrodynamical equations. 

Employing this renewed package we have created sets of hy- 
drodynamical models with different values of CR-related param- 
eters. We compare the results of our simulations with the obser- 
vations of Tycho's and SN1006 supemova remnants. 

The paper is structured as follows. In Section 2 we describe 
the basic equations we use for the simulations, in Section 3 we 
show the results of modeling. We compare our models with 
the observations in the Section 4. We summarize the results in 
Section 5, discuss them in Section 6 and conclude by Section 7. 

2. Basic equations and methiod 

2.1. Code description 

For modeling the evolution of supemova remnants, we employ 
the hydrodynamical code supremna, which was introduced by 
Sorokina et al. (2004). The method accounts for electron ther- 
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mal conduction and includes self-consistent calculations of time- 
dependent ionization processes. The electron and ion tempera- 
ture equilibration processes are parametrized. 

The code uses an implicit Lagrangean formulation for 
one dimensional spherical-symmetrical geometry. The hydro- 
dynamical evolution of the remnant is coupled with the sys- 
tem of kinetic equations of ionization balance to calculate self- 
consistently non-equillibrium ionization state of the shocked 
plasma. Ion and electron temperatures are treated separately tak- 
ing into account electron thermal conduction (see Appendix A). 

To describe the effects of collisionless energy exchange, 
Sorokina et al. (2004) introduced a parameter qt (0 < qt < 1) 
which specifies a fraction of artificial viscosity Q, added to the 
pressure of ions in the equations (A. 3), (A. 5) and plays a role 
of a source term (the details are presented in the Appendix). 
If only the coUisional exchange is taken into account, then 
qi - (1.0 - nig /nip) and the standard system of equations with 
the heating of just ions at the front is being solved. 

The artificial viscosity (Richtmyer & Morton 1967) term is 
defined as follows 



+ Anp—[r FcB.]-Anp— ?cr2, 



Q 



AqpiAuf if (Am) < 
otherwise, 



(1) 



■ velocity difi'er- 



with dimensionless parameter = 2 and Am ■ 
ence at neighboring mesh points. 

2.2. Cosmic-ray diffusion equation 

In order to tailor the scheme for simulations of SNRs we need 
to take into account cosmic-ray (CR) diffusion. We used two- 
fluid approximation (e.g. Kang & Jones 1990, Ko 1995, Malkov 
& O'C Drury 2001, Blasi 2002, 2004, Wagner et al. 2006, 
Zirakashvili & Aharonian 2010) and introduced additional CR 
diffusion equation into the (A.1)-(A.5) set. 

The one-dimensional CR diffusion equation in the Eulerian 
frame for the plane-parallel case reads 



dEcR ^ d{uEcvd 



dt 



dr 
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where Pqk = (Tcr - 1)£cr is CR pressure, iscR — CR energy 
density, is an external source of CR energy injection, /ccr — 
diffusion coefficient. We assume equation of state for CR matter 
with the fixed adiabatic exponent of ycR =4/3, 

We generate relativistic particles pressure Pcr from the arti- 
ficial viscosity Q term by introducing a parameter ^cr> that reg- 
ulates the injection of CR particles. Thus we define the source 
term as = qcRQ{du/dr) (see Appendix B). 

Note that a somewhat similar method was used by Zank et al. 
(1993), where the introduction of a source term into CR dif- 
fusion equation was performed via a "thermal leakage" mech- 
anism. Particle distribution function was divided in two parts, 
those particles that have momentum higher than a certain value 
po were treated as CR which propagate according to the CR dif- 
fusion equation. Thermal particles are energized due to the adi- 
abatic compression or anomalous heating within a subshock. 

In our study we are not concerned with microphysics of the 
generation and escape of the energetic particles, but more with 
hydrodynamical consequences of the acceleration. 

After transformation of the equation (2) to our adopted 
Lagrangian frame, we get 
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where DEcr/DI = dEcR/dt + Anr^upidEcnldm). 

We treat cosmic-ray flux Fqr in the same manner as it is 
done for thermal electron conduction in Sorokina et al. (2004). 
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where FsatcR = cPcr/2 (c — speed of light) is saturated value 
of the cosmic-ray flux (the Eddington approximation). At this 
stage we assume a constant diffusion coefficient kcr, as we do 
not have information on the spectrum of cosmic rays. 

In reality the diffusion will depend on particle energy. In that 
case the diffusion coefficient used by us should be regarded as 
an pressure weighted mean value. Moreover, according to non- 
linear cosmic ray acceleration theory, efficient acceleration leads 
to hard spectra, in which case most of the energy, and hence 
pressure, is contained by the particles with the highest energies. 
If the maximum energy is around lO'^ eV we expect a typical 
diffusion coefficient, under the assumption of Bohm-diffusion, 
of /CcR = iPrgC = 3 X 1026(B/100//G)-H£/10'^ eV) cm^s"' (r^ 
— gyroradius, B — magnetic field, E — energy of the particle). 
The largest role of the diffusion coefficient for our calculations 
is in the role of CR escape, as it drains energy from the plasma. 

And finally, we alter the equation (A. 3) by adding a compo- 
nent of relativistic particle pressure Pcr in such a way that 



du 
dt 
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dm 
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The entire system of equations (A.l), (5), (A.4), (A.5) and 
(3) describes the system, where electrons, ions and cosmic-ray 
components are treated independently. The contributions to the 
pressure of these three components are governed by two free 
parameters qi and ^cr- 

3. Numerical models 

We created a library of numerical models with various sets of 
parameters qcR and kqr. In the simulations we used a delayed- 
detonation thermonuclear explosion model with E - 1.4 x 
10^' erg (Woosley et al. 2007), initially the CR pressure is 
^CR ~ 10"^''Pi, and the temperature of the homogeneous ambi- 
ent medium is To = 10"* K. Using the artificial viscosity source 
term, we turn on CR generation only at the forward shock. 

We considered two sets of models of supernova remnants. 
Parameters for one set (Tycho's case) were taken similar to 
Tycho's SNR: the remnant is surrounded with homogeneous cir- 
cumstellar matter of density po = lO"^'* gem""* and the age of 
the system is 440 years. Note, that there are indications, that the 
real ambient density in Tycho vicinity is lower. For example, 
Katsuda et al. (2010) from proper motion measurements found 
that TiQ < 0.2 cm"-', thus po ^ 0.4 x 10"^"* g cm"-'. Nevertheless, 
taking into account that we consider over-energetic initial ex- 
plosion and that radius of the remnant R oc (E/po)^^^ is a weak 
function of ambient density, we assume that our input parame- 
ters match Tycho's SNR. 

Examples of hydrodynamical profiles of a few of these mod- 
els are presented in Fig. 1. Top row shows the simulation with 
qcR = 0.0, middle row — ^cr = 0.7, kcr = lO^^cm^s"^; bot- 
tom row — ^cr = 0.99, a-cr = 10^^ cm^s"'. Left panels: density 
(black solid, scale at the left-hand side), ion pressure (blue dash- 
dotted, scale at the right-hand side), CR pressure (blue dashed. 
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Fig. 1. HD profiles for the SNR models at f = 440 years with pcsM = lO"^'* gcm"^. Top row shows the simulation with qcR = 0.0, middle row — 
qcR = 0.7, KcR = 10"^ cm-s"' ; bottom row — ^cr = 0.99, kqr = 10^^ cm^s"' . Left panels: density (black solid line, scale at the left-hand side), ion 
pressure (blue dash-dotted line, scale at the right-hand side), CR pressure (blue dashed line, scale at the right-hand side). Right panels: velocity 
profile (black solid line, scale at the left-hand side), electron temperature profiles (blue dashed line, scale at the right-hand side), ion temperature 
profiles (dashed-dotted line, scale at the right-hand side). 



scale at the right-hand side). Right panels: soUd line shows ve- 
locity profile (scale at the left-hand side), dashed line — electron 
temperature profiles, dashed-dotted line — ion temperature pro- 
files (scale at the right-hand side). 

Physical parameters of another set (SN1006 case) were cho- 
sen to match the remnant of SN1006 conditions: an ambient ho- 
mogeneous circumstellar matter of density po = 4x10"^^ g cm"-' 
and the age of 1000 years. Hydrodynamical profiles of a few of 
these models are presented in Fig. 2. 

The numerical models were verified to satisfy the Rankine- 
Hugoniot condition, derived for the system with energy losses 



(e.g. Malkov & O'C Drury 2001, Bykov et al. 2008, Vink et al. 
2010). Namely the following relation was checked 



r- ^1 + 2(r2 - DGesc/Ooov^) 

where X - Pmax/Po is the compression ratio behind the shock, y 
— adiabatic index, Qesc — energy lost by the system, vs — blast 
wave speed. We rewrite the energy losses term from the energy 
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. Top row shows the simulation with ^cr = 0.0, middle row 
— 9cR = 0.7, KcR = 10'^ cm's"' ; bottom row — qcR = 0.99, kcr = lO"^*" cm'^s"' . Left panels: density (black solid line, scale at the left-hand side), 
ion pressure (blue dash-dotted line, scale at the right-hand side), CR pressure (blue dashed line, scale at the right-hand side). Right panels: velocity 
profile (black solid line, scale at the left-hand side), electron temperature profiles (blue dashed line, scale at the right-hand side), ion temperature 
profiles (dashed-dotted line, scale at the right-hand side). 



flux conservation law (e.g. formula (12) in Vink et al. 2010) via 
a dimensionless parameter as 

, I E + P 

fesc = 2ee,sc/(P0Vs) = 1 " ^ TTT (7) 

where E,P,p are total energy density, pressure and matter den- 
sity behind the shock. 

The top panel of Fig. 3 shows the relation (6) for 7 = 5/3 
— solid line and y = 4/3 — dashed line. Data from the nu- 
merical models are presented with colored circles. Hue reflects 
the efliciency of CR acceleration, so that red point correspond 



to the models with ^cr = 0.1 and cyan — to the models with 
^CR = 0.99. 

The bottom panel of Fig. 3 shows evolution of the compres- 
sion ratio and ratio of the contact discontinuity (CD) to blast 
wave (BW) radii. The profiles are created from the Tycho's case 
models with po = 10"^"* gcm^^*. Dashed lines show the simula- 
tion with qcR - 0.7, solid lines — ^cr = 0.9. Different colors 
correspond to diff'erent values of kqr (see legend). The curves are 
compared to those from Volk et al. (2008). We see that x pro- 
files differ from the curves presented for Tycho's SNR in Volk 
et al. (Fig. 1 in 2008). In our approach the compression ratio is 
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an increasing function of time, whereas in VoLk et al. (2008) x 
decreases with time. Radii ratio in our cases also does not follow 
exactly the same trend as in Volk et al. (2008). 

This divergence in the results can be explained by the mag- 
netic field behind the shock. Volk et al. (2005) show that on 
average the downstream magnetic field pressure B^/{8n) is pro- 
portional to preshock gas ram pressure poVg. The ram pressure 
and the magnetic field dilutes as the remnant expands, thus the 
acceleration decreases. In our models, the acceleration efficiency 
depends only on the ram pressure (Eq. 1) and does not take into 
account the vanishing magnetic field. Thus with a constant ^cR 
parameter, we somewhat underestimate the efiiciency of the CR 
acceleration at the earlier times of the evolution and we overes- 
timate it at the later times. 

Note, that for the simulation with eflicient acceleration ^cr = 
0.9 and with difliision coeflftcient kqr = 10^' cm^s"' (soUd ma- 
genta line) at the age of more than 1000 years a density spike 
starts to form. Structure of such a spike is illustrated in the bot- 
tom rows of Fig. 1 and Fig. 2. Similar spikes are also present in 
the CR dominated shocks in simulations of Wagner et al. (2009). 
Nevertheless, this type of structure is produced only at the ex- 
treme values of CR parameters and is unlikely to be reaUzed in 
nature. 

The CR precursor can be traced in the hydrodynamical pro- 
files, presented in Fig. 1 and Fig. 2, but the precursor region in 
front of the blast wave is thin (~ 10^^ cm) and cannot be accu- 
rately resolved. 

4. Comparison with the Galactic SNRs 

There is an evidence that the shocks of Tycho (e.g. Warren et al. 
2005) and SN1006 (Cassam-Chenai et al. 2008, Orlando et al. 
2008, Miceli et al. 2009) are modified by acceleration and diffu- 
sion of CR particles. Thus, we applied our numerical models to 
these Galactic supernova remnants. Locations of the blast wave 
(BW), the contact discontinuity (CD) and the reverse shock (RS) 
measured by Warren et al. (2005) for Tycho's SNR, were com- 
pared with the corresponding values obtained in our simulations. 

The right plot of Fig. 4 shows the measured (horizontal strip) 
radii ratios and also results of the simulations (asterisks). To ac- 
count for the three-dimentional projection effects and Rayleigh- 
Taylor (R-T) instabiUties^ that can produce fingers of ejecta that 
protrude out toward the BW (Chevalier et al. 1992, Wang & 
Chevaher 2001) we adopt a correction factor of 1.07 (Warren 
et al. 2005, Cassam-Chenai et al. 2008, Volk et al. 2008, and ref- 
erences therein) to the modeled CD:BW radii ratios. The models 
that match the observed radii are boldfaced. 

We also performed simulations with CR acceleration at the 
reverse shock with the same injection efficiency ^cr and diffu- 
sion coeflScient ^cr as at the blast wave. In these models the layer 
of the shocked supernova ejecta becomes thinner, and the reverse 
shock approaches too close to the forward shock compared to the 
data from Tycho's SNR. Thus, we did not study these scenarios 
in details. 

From an analysis of the data, presented in Fig. 4, we derive 
that the most plausible models for Tycho are with qcR - 0.5 - 
0.7, KcR = lO^"* - 10^5 cm^s-'. In fact, models with ^cr = 0.3, 
""CR = 10^*" cm^s"' and ^cr = 0.9, kqr = 10^^ cm^s"' also 
match the observation. Nevertheless, if we assume that diffusion 
coefficient kcr increases with the energy of the CR particles and 
that the amount of energetic particles increases with the injection 

' Note, thai Ihc R-T instability is not considerably affected by accel- 
eration at the forward shock (Ferrand et al. 2010, Wang 2010) 



efficiency qcR, then we conclude that low efficiency and high 
diffusion or high efficiency and low diffusion scenarios probably 
are not realized. 

For Tycho, the simulations yield the compression ratio x = 
4.3 - 6.8. These results are in agreement with Cassam-Chenai 
et al. (2007) and Volk et al. (2008). 

The similar comparison of the models to with SN1006 rem- 
nant data (Miceli et al. 2009) is presented on the right plot of Fig. 

4. From these models we obtain the compression ratio behind 
the blast wave of the remnant = 4.7 - 8.3. All the models with 
^CR = (0.3-0.9) match the observational data. Note, that models 
of SN1006 with constant adiabatic index y (Petruk et al. 2011) 
could not explain the observed small distance between forward 
shock and contact discontinuity. 

5. Results 

We developed a hydro-code to simulate an evolution of a 
spherically-symmetrical supernova remnants with account for 
CR acceleration and time-dependent ionization. We created two 
sets of hydrodynamical models with different values of pa- 
rameters for CR acceleration efficiency and CR energy losses. 
The comparison of these models with the measurement of 
BW:CD:RS radii of the Galactic SNRs suggests the following. 

The simulations with ^cr = (0.4 - 0.7) and kqr = (10^"* - 
lO^^") cm^s ' match the observed radii ratios of Tycho's SNR 
(Fig. 4). This range of values covers the estimate kcr = 2 x 
lO^'* cm^s"' found by Wagner et al. (2009) for Tycho and are in 
agreement with the findings of Parizot et al. (2006) and Eriksen 
etal. (2011). We derived the compression ratio of = (4.3-6.8), 
energy losses due to CR difftision is of Cesc = (0.1 - 0.2) and 
distance to Tycho's SNR of 3.3 pc ^. The distance is in agreement 
with the results of Hayato et al. (2010) and Tian & Leahy (201 1). 

For SN1006 remnant we obtained = (4.7 - 8.3) with losses 
of fesc = (0.2 - 0.5). We derived the distance to the SN1006 
remnant as 2.0 kpc. 

6. Discussion 

The observed ratios of SNR reverse shock, contact discontinuity 
and forward shock radii may depend on a number of yet un- 
known factors, such as structure of the CSM (potential presuper- 
nova wind and density enhancement, see discussion in Kosenko 
et al. 2010, Xu et al. 2010), CR acceleration and diffusion effi- 
ciency. Thus the results, presented in this study do not claim to 
be exhaustive, but rather describe one of the possible scenarios 
which explains the observations within the cosmic-ray accelera- 
tion paradigm. 

In the numerical models considered here, we do not tum on 
CR acceleration at the reverse shock. Results from the simu- 
lations with the same acceleration efficiency (<7cr) of the rel- 
ativistic particles at the reverse shock do not match the ob- 
served Rrs/Rbw ratio for Tycho's SNR. Helder & Vink (2008), 
Zirakashvili & Aharonian (2010) point out to the possible ac- 
celeration of particles by the reverse shock of CasA and SNR 
RX J1713. 7-3946. In addition, numerical simulations of Schure 
et al. (2010) showed, that the reverse shock can re-accelerate 
particles, provided that the magnetic field is sufficiently strong. 
However, our models indicate that for Tycho's SNR the cosmic 



^ This value depends on the assumed explosion energy and CSM den- 
sity (e.g. decrease of po by a factor of three results in increase of the 
distance estimate by 25%). 
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Fig. 3. Top panel: Rankine-Hugoniot relation (6) and data points from all the numerical models. Hue reflects the efficiency of CR acceleration. 
Bottom row: left panel shows evolution of compression ration ^, right panel shows evolution of contact discontinuity to blast wave radii ratio, 
Po = 10"^'* g cm"^. Dashed lines show models with ^cr = 0.7, solid lines — ^cr = 0.9. Different colors indicate different values of the diffusion 
coefficient (see legend, where kqk values are given in units of cm^s"'). 



rays acceleration at the reverse shock is not as efficient as at the 
forward shock. 

The possible factors which suppress the acceleration at the 
reverse shock are as follows. The reverse shock velocity in the 
frame of ejecta is lower in comparison with the blast wave, and 
acceleration efficiency is probably an increasing function of the 
shock velocity. Moreover, the initial magnetic field of the pro- 
genitor white dwarf is diluted in the remnant by many orders of 
magnitude. Non-vanishing magnetic field is necessary for the ac- 
celeration to take place. Also, low efficiency of CR acceleration 
at the reverse shock was found for Kepler SNR by Decourchelle 
et al. (2000). 

We compared the input and output parameters of the simula- 
tion with Rankine-Hugoniot (RH) relations, derived in two-fluid 
approach for the steady-state situation and a plane-parralel ge- 
ometry, presented in Vink et al. (2010). On the one hand, the 
values of PcR/^tot measured in the models directly behind the 
shock cannot be compared to w defined in Vink et al. (2010). 
The diffusion coefficient a-cr implements the CR energy losses 
and thus the CR pressure Pcr is a non-monotonic function of /ccr 
(Fig. 5). Namely, the pressure decreases with kcs. in the regime 
of efficient acceleration, where a-cr > 3 x 10^^ cm^s"'. This 



is opposite to the behavior of w defined by Vink et al. (2010). 
Possible explanation is that spherical symmetrical expansion and 
inward CR diffusion efficiently dilute the pressure, while in the 
plane-parallel case described by Vink et al. (2010) these effects 
are absent. 

On the other hand, the parameter ^cr sets the share of the 
CR pressure taken from the entropy pool created by the shock 
and in some sense can be attributed to w. Fig. 6 shows profiles 
of compression ration eesc (left panel) and CR energy losses x 
(right panel) as a function of w = P cr IP tot (Fig. 1 in Vink et al. 
2010, for Mach number Mo = 500). In the same panels we plot 
eesc and X measured in the models versus the input parameter 
^CR- The intensity of the data points corresponds to the value 
of a-cr: the lowest value corresponds to the black asterisks, the 
models with the highest value of a-cr shaded in light gray. 

The right-hand plot of Fig. 6 can be explained with the help 
of the light-hand plot of Fig. 5, where compression ratios versus 
diffusion coefficient for different acceleration efficiencies qcR 
are presented. This plot shows that x has a maximum at cer- 
tain A-cR. These maximal values of x for different ^cr follow the 
RH curve, plotted with thick black line in the right plot of Fig. 6. 
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9cR (Tycho's case). 



Respectively, the maximum attainable energy losses are skirted 
by the RH curves, presented in the left panel of Fig. 6. 

The parameters ^cr and kcr are fixed for each model in our 
method, while in reality they vary with time and location (e.g. 
Zirakashvili & Aharonian 2010). The effects of the variability 
of the diffusion coefficient was also discussed extensively by 
Lagage & Cesarsky (1983), where they point out that kqr in- 
creases away from the shock. Nevertheless, we assume that it is 
acceptable to use averaged constant values in this approach. 

Even though we used two parameters for CR component in 
our models, in fact, the theory of diffusive CR acceleration as- 
sumes that, <7cR and kcr are not independent. Thus it is possible 
to incorporate a physical relation between those two in future, 
reducing the problem to a case with only one free parameter. 

At last we note, that in comparison of the models and the ob- 
servations, we considered also a case with lower R-T correction 
factor of 1.05. In this case the models that match the observed 
radii yield the compression ratio for Tycho's SNR systematically 
higher with;t' - 4.6 - 8.5, but within the uncertainties of our ap- 



proach and errors of radii measurements, the final results are not 
altered considerably. 



7. Conclusion 

We studied hydrodynamical models of supernova remnants evo- 
lution with account of diffusive cosmic -ray acceleration. We ap- 
plied a two-fluid approach to investigate the effects of acceler- 
ation efficiency and energy losses on observable properties on 
the remnant. We compared the results of our simulations with 
the measured radii of BW, CD and RS of Tycho's and SN1006 
SNRs. 

We analyzed the numerical models and checked that the 
Rankine-Hugoniot relation for compression ratio and energy 
losses is met (Fig. 3). This method has a list of shortcomings but 
it includes all important relevant physical processes and easy to 
use for fast modeling of supernova remnants and for comparison 
with observations. 
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Fig. 6. Left panel: escape CR energy versus w and t/cR- Right panel: compression ratio ;t- versus w and ^cr (Tycho's case). The Rankine-Hugoniot 
black curve corresponds for the case of Mach number Mq = 500. Intensity of the data point is invers proportional to the values of diffusion 
coefficient /Ccr, i.e. black asterisk correspond to Kcb. = 10^ cm^s"', light grey correspond to Kqr = 10^' cm^s"^ 



We found that, in order to explain radial properties of 
Tycho's SNR, the simulations yield the compression ratio be- 
hind the blast wave of (4.3 - 6.8), energy losses due to cosmic 
ray particle diffusion of (0.1 - 0.2)poVg/2. Distance to Tycho is 
estimated as 3.3 pc. For the case of SN1006 remnant, we found 
X = (4.7 - 8.3) with losses of (0.2 - 0.5) pov^/2. The distance to 
the remnant is 2.0 kpc. 

In the future we plan to employ the developed package for 
calculation of detailed thermal X-ray emission from the super- 
nova remnant models modified by the acceleration. The resulting 
synthetic spectra will allows us to study the effects of different 
physical conditions of the shock plasma on the X-ray spectra. 
The analysis will be performed for different explosion models 
(similar to the method of Badenes et al. 2006, 2008) and later 
compared with observations of young remnants of type la super- 
nova. 
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Appendix A: Basic equations 

Originally the method (described by Sorokina et al. 2004) solves 
the following system of hydrodynamical equations 
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(A.6) 



Here u is the velocity; p is the density; Tg and T, are the electron 
and ion temperatures; Pe and f , are the respective pressures, tak- 
ing into account artificial viscosity (see description below); Eg 
and Ej are the thermal energies per unit mass of the gas element 
at the Lagrangean coordinate m (corresponding to mass within a 
radius r) at the moment of time t; F^ond is the energy flux due to 
the electron-electron and electron-ion thermal conduction; v,g is 
the electron-ion collision frequency per unit volume; Sr is the ra- 
diation energy loss rate per unit mass of the gas element; dsion/dt 
accounts for the change in specific thermal energy of gas due to 
the change of ionization state; X = {Xhi, Xun, Xuei, ■ ■, ^Nixxix} 
is the abundance vector of all ions of all included elements rela- 
tive to the total number of atoms and ions. Also X^ = ng/ri/, for 
the number of electrons per baryon is introduced. The physical 
processes described by these equations are discussed in Sorokina 
et al. (2004). 

Appendix B: On the artificial viscosity 

Originally in supremna the ion pressure was generated as P, = 
f,(phys) -I- qiQ, thus for the electron pressure we put Pg = 
/'e(phys) -I- (1 - qi)Q. Sign "(phys)" stands for physical pressure 
in the system. 

The meaning of this pressure can be clarified from the fol- 
lowing simple thermodynamical relation. The second law reads 
as dE+PdV = 0, with E — intemal energy, P — pressure, dV — 
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volume element. Rewriting the pressure term as f = f (phys)+2, 
we get dE + (P(phys) + Q)dV = 0. Thus, the term Q plays a role 
of a source function in dE + P(phys)<iV = -QdV. If there are 
energy losses e in the system, then the final relation reads as 
dE + PdV = -QdV + s. 

After the introduction of the CR component, using ^cr pa- 
rameter, we define Pgr = PcR(phys) + qcviQ (^cr = 0.0 de- 
scribes a case where the contribution from the cosmic-ray com- 
ponent is zero). Thus, we have the following distribution of the 
artificial viscosity (Eqn. 1) between ion and electron pressure: 
Pi = P;(phys) + (1 - qcsdmQ and = P.Cphys) + (1 - - 
qi)Q. Therefore equation (2) yields the source term in the form 
oi@^qcRQ{duldr). 
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